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Abstract. - Reverse engineering of complex dynamical networks is important for a variety of fields 
where uncovering the full topology of unknown networks and estimating parameters characterizing 
the network structure and dynamical processes are of interest. We consider complex oscillator 
networks with time-delayed interactions in a noisy environment, and develop an effective method 
to infer the full topology of the network and evaluate the amount of time delay based solely 
on noise- contaminated time series. In particular, we develop an analytic theory establishing 
that the dynamical correlation matrix, which can be constructed purely from time series, can be 
manipulated to yield both the network topology and the amount of time delay simultaneously. 
Extensive numerical support is provided to validate the method. While our method provides a 
viable solution to the network inverse problem, significant difficulties, limitations, and challenges 
still remain, and these are discussed thoroughly. 
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Time-delayed interactions are common in com- 
plex systems arising from various fields of science 
and engineering. Con- sider, for example, a cou- 
pled oscillator network in a physi- cal environment 
where noise is present. Time delay can typically 
occur in the node-to-node interactions. Now sup- 
pose that no prior knowledge about the nodal dy- 
namics and the network topology is available but 
only a set of noise- contaminated time series can be 
obtained through measure- ments. The question 
is whether it is possible to deduce the full topol- 
ogy of the network and to estimate the amount 
of average time delay using the time series only. 
This issue belongs to the recently emerged sub- 
field of research in com- plex systems: reverse 
engineering of complex networks (or the inverse 
problem). While a number of methods address 
the network inverse problem have appeared, to our 
knowl- edge, the issue of time-delayed interactions 
has not been considered. Here we present an effec- 
tive method to infer the full network topology and, 



at the same time, to estimate the amount of aver- 
age time delay in the network. In partic- ular, we 
develop a physical theory to obtain a formula relat- 
ing the network topology and time delay to the 
dynamical correlation matrix, which can be con- 
structed purely from time series. We then show 
how information about the time delay encrypted 
in the dynamical correlation matrix can be sepa- 
rated from that of network topology, allowing both 
to be inferred in a computationally extremely ef- 
ficient man- ner. We present numerical examples 
from both model and real-world complex networks 
to demonstrate the working of our method. Dif- 
ficulties, limitations, and challenges are also dis- 
cussed. Reverse-engineering of complex dynamical 
systems has potential applications in many disci- 
plines, and our work represents a small step for- 
ward in this extremely challenging area. 

One of the outstanding issues in nonlinear and statisti- 
cal physics, and also in network science and engineering, is 
to infer or predict the topology and other basic character- 
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istics of complex networks based only on measured time 
series. This "inverse" problem is relevant to a number of 
fields such as biomedical and techno-social sciences where 
complex networked systems are ubiquitous. In defense, 
the problem of identifying various adversarial networks 
based on observations is also of paramount importance. 
Despite tremendous efforts in revealing the connection be- 
tween network structures and dynamics [IHS] , how to infer 
the underlying topology from dynamical behaviors is still 
challenging as an inverse problem, especially in the ab- 
sence of the knowledge of nodal dynamics. 

Recent years have witnessed the emergence of a number 
of methods to address various aspects of the inverse prob- 
lem, which include gene networks inference using singular 
value decomposition and robust regression [6], spike clas- 
sification methods for measuring interactions among neu- 
rons from spike trains [7] , symbolically reverse engineering 
of coupled ordinary differential equations [5] , approaches 
based on response dynamics of a specific oscillators [§], L\ 
norm in optimization theory [10] . noise induced scaling 
laws [11] . and the interplay between dynamical correla- 
tion and network structure in the presence of noise [12] . 
However, the issue of time delay has not been addressed 
yet. The purpose of this article is to present a general the- 
ory that leads to a completely data-driven and extremely 
efficient method to predict the network topology and the 
time delay at the same time. 

Time delay is fundamental in natural systems, due to 
the finite propagation speed of physical signals. In ad- 
dition to numerous examples in physics, situations where 
time delay is important include the latency times of neu- 
ronal excitations in neuroscience, finite reaction times of 
chemicals in chemistry, etc.. In coupled oscillator net- 
works, the effects of time delay on dynamics under vari- 
ous given network topologies have been studied extensively 
[13H17] . In our case, however, the network connections, 
the amount of the time delay, and other properties of the 
network are unknown a priori, and our goal is to predict 
these by using noisy time series only. To be concrete, 
we shall focus on complex oscillator networks. Our gen- 
eral point of view is that, information about the network 
topology and time delay has been encoded in noisy time 
series from various nodes in the network. The objective of 
solving the inverse problem is to decode such information 
from noisy time series. 

Our idea is that, if the networked system suffers from 
a noisy environment so that the measured time series are 
noisy, it is possible to accomplish the task of decoding in a 
natural way. In particular, we construct a dynamical cor- 
relation matrix from all available time series, the elements 
of which are the average products of the deviations of all 
pair-wise time series from a mean value. We shall show 
analytically that information about the network structure 
and time delay can be decoded through this matrix. In 
fact, as we will show in developing our theory, informa- 
tion about the network topology can be separated from 
the time delay through a generalized inverse operation of 



the dynamical correlation matrix, enabling a complete pre- 
diction of the underlying networked system. 

To provide numerical support for our theory, we exploit 
three representative dynamical systems on homogeneous 
and heterogeneous model complex networks and on a num- 
ber of real-world networks as well. We find that the pres- 
ence of a time delay results in a deviation in the distribu- 
tion of the diagonal elements of the dynamical correlation 
matrix from a power law, which can be used as a prelim- 
inary criterion to determine whether there is a significant 
time delay in the underlying networked system. Compu- 
tations reveal high accuracies in the prediction of both the 
network topology and the time delay for all combinations 
of dynamical systems and network models studied. 

We present our theory and method by considering a 
network of ./V coupled oscillators. Each oscillator, when 
decoupled, satisfies x; = Fj[xJ, where Xi denotes the d- 
dimcnsional state variable of node i. The dynamics of 
the whole time-delayed system in a noisy environment is 
described as: 

N 

ki(t) = Fj[xi(t)] - c^LyH^fi - r)] + ffi(t), (1) 

where c is the coupling strength and H denotes the cou- 
pling function. Lij is Laplacian matrix, characterizing the 
topology of the underlying network that L,j = —1 if j 
connects to i (otherwise 0) for i ^ j, and La = ki, where 
ki is the degree of node i. r denotes the time delay, and ffi 
is a d-dimcnsional stochastic process representing noise on 
node i (In the following, we use ~* on the head to denote 
the <i-dimensional state variable). The standard proce- 
dure of linearization PQQ2] can be carried out by letting 
Xj = 5Q+£i, where x^ is the counterpart of x, in the absence 
of noise. The d-dimcnsional dynamical process governing 
the fluctuations on ith oscillator can then be obtained as 
the variational equation: 

ii (t) = DFj ■ & (t)-cJ2 Lij DU ■ £• (t-T)+fk(t), (2) 
i=i 

where DFi and DU denote the d x d Jacobian matrices 
of the intrinsic dynamics Fj and the coupling function H, 
respectively. Decomposing Eq. ([2]) in terms of the eigen- 
modes, we obtain 

= Y, D¥ "P ■ e ~M*) _ cA «^ H • - T ) + £»(*)■ (3) 

p 

Here, instead of the index i,j running on the real space 
of networks, the index a,/3 run on the eigen-space. e a = 

Y^i^aid, (a = 52il/>aifk and D¥ af} = J^i^aiDFilppi, 

where tp a j denotes the ath normalized eigenvector of the 
Laplacian matrix, and A Q is the corresponding eigenval- 
ues that satisfy = Ao < Ai < • • • < Aiv-i- Under the 
approximation DFi ~ DF so that -DF Q( g = DFSap, the 
above equation can be reduced to 

r„(t) = DF ■ e a (t) - c\ a DH ■ e a (t - r) + £*(*)• (4) 
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a 2 IdSijS(t — t') with Id the d-dimensional identity ma- 



From the covariance of Gaussian noise (rfi(t)f)j (t 1 )) = time delay r and the structure information in terms of 

= J2 a =i ^aiipaj/^a, the pseudo- inverse of Laplacian 
matrix. Note that the time delay has no effect on the cross- 
correlation elements except the auto-correlations due to 
the identity matrix. Following Eq. ([5]), the diagonal ele- 
ments Cu of the dynamical correlation matrix can be ob- 
tained by expanding in terms of the underlying network 
structure [12] : 



trix and a the noise strength, we obtain (£ Q (i)£i (t')) = 
<r 2 IdS a i3 8(t~t'), which indicates the stochastic process we 
mapped into eign-space is still Gaussian noise. Assuming 
small time delay, we can apply the first-order approxima- 
tion: e a (t — r) = e a (<) — re a (t) , which yields 



(I d - cT\ a DK)e a (t) = -{c\ a DH - DF)e a {t) + £„(*). 

Denote B = (I d - ctA q £>H) _1 , A = B(cX a DH - DF), 
and follow the standard stochastic calculus [TS], we get 
the solution: 



Ea(*) = 



-A/ 



I e-^'-^BCaW- 
Jo 



(5) 



Since we are interested in the regime where oscillator 
states are perturbed from the synchronized manifold by 
the noisy environment, we assume the system is in the ab- 
sence of divergence of state variables, therefore in the long 
time limit, the initial condition term can be discarded and 
we have |12j : 

A(e a el) + (e a el)A = er 2 BB T . 

The general solution of {e a <^), the d x d covariance ma- 
trix about ci-dimensional states of the ceth oscillator in the 
eigen-space, can be written as [19j : 

vec((e a el)) = cr W(BB T )/(I d ® A + A <g> I d ), 

where the operator uec(X) creates a column vector from a 
matrix X by stacking the columns of X below one another. 

Although we obtain this solution, it is not practical in 
real applications. In follows, we approximate the state 
variables as one dimension such that DH = 1 and drop 
the notation^. In this way, the above solution is simplified 
to: 



(4) 



1 



2c (1 - CT\ a )(\ a - DF/c)' 



(6) 



Return to real variables from the eigen-space by inserting 
Ci = Ea ^«! £ « m t° the correlation function CV, = 
in the real space between any two nodes, we have C\j = 



J2a=l V'aiV'aj^a) SUCn tnat 
9 N-l 



Ci 



-T- 

2c ^ (1 



J (l-cr\ a )(X a -DF/c) 



(7) 



Under the approximation of negligible DF/c and remind- 
ing of the small time delay r, Eq. ([7]) for the dynamical 
correlation can then be expanded as: 

2 N - 1 1 , \ 2 

^ a \ ^ 1 + cr\ a a rT + . , . 

Ci i ~ 2c~ 1^ \ Vaii' a] = — [L 1 + crljvjij, (8) 

a—l a 

wherein under the effect of noise a 2 , the dynamics in terms 
of correlation matrix C is connected explicitly with the 
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2c 
2 



[K 1 + K : PK 1 + K^PK^PK" 1 ] 

a / 1 



a- 



1+ m 



where K = diag(fci, ■ ■ ■ , fcjv) is the degree matrix, P is the 
adjacency matrix such that L = K — P and (k) denotes 
the average degree. We see that the fluctuations Cu at 
node i depend both on its local structure hi and the time 
delay r. When r = 0, this result is consistent with the 
recently discovered noise- induced scaling law [TT] , derived 
there by a power-spectral analysis. 

The off-diagonal elements of L contain complete infor- 
mation about the network structure while its diagonal el- 
ements can be obtained from the off-diagonal ones. We 
thus focus on the off-diagonal elements. For i ^ j, follow- 
ing Eq. ([7]), the generalized inverse matrix is 



CI 



2c 



N-l 



r 2 



\ a {l-CT\ a )ll) a i1p 



2c 



[L-ctL 2 ] v , (10) 



Considering L = K — P, we can cast this equation in the 
following form: 



CI 



2c 



2 [L + cr(KP + PK K 2 P 2 )],, . 



(11) 



For those off-diagonal elements the diagonal matrix 

K 2 has no contributions and P 2 contributes lij, where 
Uj is the number of two-step paths connecting i with j. 
By considering the negligible contribution of lij compared 
with degrees, we thus have 



2c ^ 



Lij + cr(ki + kj), if i connects with j 



0. 



otherwise 



(12) 



Equation (|12[) is one of our main results for network in- 
ference, which indicates that the network structure can be 
inferred through the off-diagonal elements Cjj of the dy- 
namical correlation matrix based solely on the measured 
time series. 

Once L is predicted, the time delay r can be estimated, 
e.g., from Eq. (fTU|) . We obtain 



oMij 



(13) 



i^,i» 3 ^0,(L 2 ) i3 ^0 



where the subscript in the average (•) covers all possible 
pairs of i and j by excluding the diagonal elements in the 
matrices L and L 2 , and all pairs with zero elements in the 
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Fig. 1: (Color online) Diagonal elements Cu of the dynami- 
cal correlation matrix as a function of node degree k for three 
dynamical processes with different time delay r on scale-free 
and random networks. Square, circle, triangle and reverse tri- 
angle denote r = 0.01, 0.05, 0.07 and 0.09, respectively. The 
curves are the theoretical prediction from Eq. |[9j). The sizes 
of model networks are 100 and the average degree is 10. The 
noise strength a 2 is 0.1 and the coupling strength c is 0.2. 



matrix L or L 2 . Excluding zero elements can effectively 
reduce the estimation error for r. 

We now demonstrate numerically our method by con- 
sidering several model and real- world networks in the pres- 
ence of noise and time delay. For each network, we imple- 
ment three dynamical processes: (i) Consensus dynamics 



N 

i(t) = c^Pij[xj(t - t) -Xi{t-r)] +T}i 

3 = 1 



(ii) Rossler dynamics [2Tj : 
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±i(t) = ~yi - Zi + cY^ = i Pij[xj(t - t) - xt(t - t)] 
ih = x,i + 0.2y. t + cY, j=1 Pij(Vj ~ Vi), 
Zi = 0.2 + Zi { Xi - 9.0) + cY% =1 Pijfa ~ 

and (iii) Kuramoto phase oscillators [2"2"] : 

JV 

8{t) = Ui + c^sin[^(i -r) - 9 t (t - r)] + Vt , 

where 9i and cjj are the phase and the natural frequency 
of oscillator i. 

Time series are then collected from all nodes. The el- 
ement of the dynamical correlation matrix between two 
arbitrary nodes i and j is calculated as CV,- = ([xi(t) — 
x(t)} ■ [xj(jb) - x(t)]) t , where x(t) = (l/A0£*Li*i(*) and 



Fig. 2: (Color online) (a) Example of the distribution of the 
values of elements of the generalized inverse C 1 ^ of the dynam- 
ical correlation matrix C for consensus dynamics associated 
with a scale-free network, where r = 0.05. The success rate 
of prediction of existent links S e for (b) consensus dynamics, 
(c) Kuramoto oscillators and (d) Rossler dynamics as a func- 
tion of time delay r for a number of model and real-world 
networks: scale- free networks (scale- free) [23], random net- 
work (random) [24] . small- world network (small- world) [25] . 
dolphin social network (dolphins) [26], network of American 
football games among colleges (football) [27], friendship net- 
work of karate club (karate) 28 and network of political book 
purchases (book) [29] • Other parameters are the same as in 
Fig. [T] The success rate of nonexistent links is higher than 
0.99 for all considered cases and thus are not shown. 



(■)t denotes the long time average. For the Rossler dynam- 
ics, Xi(t) stands for the x component of the ith oscillator 
and for the Kuramoto dynamics, Xi(t) stands for the phase 
variable 6i(t) of the ith oscillator. 

Figure [1] provides an example of the dependence of fluc- 
tuations Cu on the time delay for three dynamical pro- 
cesses on both heterogeneous and homogeneous networks. 
The results are in good agreement with the theoretical 
prediction from Eq. ©. A non-ignorable deviation from 
the predicted fluctuations in noisy Rossler dynamics with 
tfcriiei delay is observed, which is mainly due to the simplifi- 
cation of one-dimensional variable to get Eq. (0) while the 
state variable in Rossler dynamics has three dimension. 
Note that, in the absence of time delay, the dependence 
of Cu on the node degree hi can be described as a power 
law [TT][T2]: Cu ~ fc" 1 , regardless of the dynamics. For 
t / 0, the deviation of Cu from the power-law behavior 
can then be used to assess preliminarily whether there is a 
significant time delay in the underlying networked system: 
a more severe deviation suggests a larger value of the time 
delay. 

After calculating the dynamical correlation matrix C, 
we can infer the details of the network connections through 
Eq. (fT2"j) via the generalized inverse of C. Figure r2Ja) 
shows the distribution of the off-diagonal elements of 
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Fig. 3: (Color online) Predicted time delay r' from Eq. (|13|) 
versus the true (pre-assumed) values for the three dynamical 
processes on a number of model and real-world networks. The 
symbols denote the same networks as in Fig. [5] The lines are 
t' = r. Other parameters are the same as in Fig. [1] 



[cr 2 /(2c)]Ct. We observe a bimodal behavior with two 
peaks: one centered at a negative value which corresponds 
to existent links, and another centered at zero which indi- 
cates non-existent links. Without time delay, the hump in 
the distribution for the existent links should be centered 
at —1. While with time delay, due to the contribution of 
the term cr(ki + kj) ~ 2cr(fc), the center of the hump 
will shift toward zero. The larger the time delay, the more 
significant the shift will be. For example, as shown in 
Fig. [Ha), the amount of the shift in the negative peak 
is 2cr(k) = 0.2. Considering c = 0.2 and (k) = 10 in 
the example, we obtain t = 0.05, which is exactly the pre- 
assumed value of the time delay in the system. To separate 
the two humps, a threshold is needed, where all existent 
links in the network are identified by the elements in 
which lie below the threshold. In particular, the subscript 
ij of a picked element below the threshold indicates a link 
between nodes i and j. We set the threshold at the si- 
tus, which corresponds to the minimum value of the fitted 
curve between two humps. The performance of our pre- 
diction method can be characterized by the success rate 
S e of the existent links, which is defined as the ratio of 
the number of successfully predicted existent links to the 
total number of existent ones. As shown in Figs. HJb)-(d), 
our method yields high success rates for different values of 
the time delay r, regardless of the nodal dynamics and of 
the network structures. 

After the network topology L is predicted, we can esti- 
mate the time delay r through Eq. ([13)) . As demonstrated 
in Fig. [3J the predicted values of r are quite close to the 
real values for almost all dynamical processes and network 
structures considered. The non-ignorable deviation occur- 
ring in Kuramoto dynamics, may be due to its sin coupling 
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Fig. 4: (Color online) (a) The success rate of prediction of 
existent links S e as a function of the average time delay r for 
different range of time delays for random consensus networks, 
(b) The predicted average time delay r' versus the original time 
delay r for different ranges of time delay. The lines are r' = r. 
Other parameters are the same as in Fig. [TJ 



function which is conditionally approximated to the linear 
coupling function used in our theory. 

We also examine the validity of our method for pre- 
dicting coupled network system with inhomogeneous time 
delay. A random consensus network associated with ran- 
dom time delays within a certain range is considered. The 
success rate S e as a function of the average time delay t 
among all pairs of nodes for different ranges of time de- 
lay is shown in Fig. |4|a). We see high success rate if the 
average time delay is not too large. Fig. @Jb) shows the 
predicted average time delay r' versus the original time 
delay r for different ranges of time delay. The predicted 
t' is in good agreement with r. These results demonstrate 
that our approach is applicable to interacting units asso- 
ciated with inhomogeneous time delay. 

While our theory and the prediction method are based 
on the system model Eq. (JlJ , a similar theory and method 
can be developed for variants of the model. For instance, 
one can consider the following system: 



A' 



(t)=F i [S i (t)]-c^P ij (U[x i (t)]-ll[x j (t-T)]) 



3=1 



Vi 



(14) 

with the Laplacian matrix in Eq. ([!} replaced by Pjj, 
the adjacency matrix of the underlying network. The dif- 
ference is that, for the dynamics of a given node in Eq. 
(|14[) . the time delay occurs only for the state information 
transmitted from its connected neighbors other than the 
dynamics of itself. Using similar analytical treatment, one 
may arrive at 

2c, 



C 



r[(cr(fc) - 1)P -ctP* 



and 



2c 



c+ 



(15) 



(16) 



\ c[(fc)P - P 2 b 7 - / 

The network structures and time delay can again be pre- 
dicted for various nodal dynamics and network structures, 
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as we verified through extensive numerical simulations 
(data are not shown here). 

In summary, we have established a theory to address 
the inverse problem for complex networked systems in the 
presence of time delay and noise, based solely on measured 
time series from the network. Especially, we have obtained 
a formula relating the generalized inverse of the dynamical 
correlation matrix, which can be computed purely from 
data, to the structural Laplacian (or adjacency) matrix 
and the amount of time delay. Under reasonable approxi- 
mations the network topology and the effect of time delay 
can be separated, leading to a computationally extremely 
efficient method for inferring the network topology and 
for estimating the time delay. The validity of the method 
has been tested numerically using a variety of combina- 
tions of nodal dynamics and network topology, including 
a number of real-world network structures. Our method 
is completely data driven, and we expect it to be applica- 
ble to the critical network inverse problems in a variety of 
fields, such as biomedical and social sciences where com- 
plex networked systems are ubiquitous. 

JR thanks Dr. Gang Yan for useful discussions. WXW 
and YCL arc supported by AFOSR under Grant No. 
FA9550-10-1-0083. 
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